A FINITE DIFFERENCE METHOD FOR PIECEWISE 
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Abstract. In this paper the numerical approximation of solutions of Liouville-Master Equations 
for time-dependent distribution functions of Piecewise Deterministic Processes with memory is con- 
sidered. These equations are linear hyperbolic PDEs with non-constant coefficients, and boundary 
conditions that depend on integrals over the interior of the integration domain. We construct a finite 
difference method of the first order, by a combination of the upwind method, for PDEs, and by a 
direct quadrature, for the boundary condition. We analyse convergence of the numerical solution for 
distribution functions evolving towards an equilibrium. Numerical results for two problems, whose 
analytical solutions are known in closed form, illustrate the theoretical finding. 
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1. Introduction. We deal with the following system of PDEs: 

(1.1) dtF,{x, y, t) + A,{x) d^F,{x, y, t) + dyFs{x, y, t) = -A,(y) F,(x, y, t) 
with Cauchy initial conditions: 

(1.2) F,{x,yM)^FoAx)5{y) 
and boundary conditions: 

S nt-to 

(1.3) F,(x,0,t) = Vg,, / dyF,{x,y,t)\,{y) 

, = 1 ^0 

for the s = {1, . . . , S} unknowns Fs -.V ^M., with {x, y,t) €V := (n x [0, T - to] x 
[to,?"]) C R^, where ft — [il.a,^b\ C M. qsj are the elements of a stochastic matrix 
having the following fundamental properties: < Qsj < 1 and Isj — 1- The known 
functions ^0,^(2:), As{x) and Xs{y) > 0, will be discussed later. 

Eq. (jl.ip . jointly with boundary conditions, is the Liouville-Master equation^ for 
the probability distribution functions Fs{x, y, t) of a continuous piecewise-deterministic 
process (PDF), that has been introduced by Davis [SHE]. Indeed, here we deal with a 
simplified version of Davis' PDPs, but still enough general to cover many interesting 
models. The definition of PDP is more popular between researchers working on op- 
erations research and probability calculus (see, e.g., [21), rather than others outside 
these fields, even though the latter unknowingly use it, at least in a simplified form. 
Before to proceed with the discussion of the numerical solution of our problem, we 
give a short introduction of the underling PDP process we are considering here. □ 

Definition 1.1. We name X{t), X ■.M.^M., be a continuous PDP if: 
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^In general, equations for density probability of random processes are derived from the Chapman- 
Kolmogorov equation. As discussed in Ref. [3] the same equation turns into a Liouville equation in 
absence of randomness, and into a Master Equation, if only jump processes are involved. We use 
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(a) X{t) satisfies the equation: 

(1.4) X{t)^AsiX), s^l,...,S 

where : M ^ M is a function chosen randomly on a set of {Ai, . . . , As} 
known functions. Given A^, we say that the dynamics is in the (deterministic) 
state s. We require that Ag{x) be Lipschitz continuous, so that, for fixed s, 
X{t) exists, is unique and non-explosive solution. 

(b) The initial condition is settled by the Cauchy problem to Eq. Jj.^l ), i.e. 
X{to) = Xo, and by the initial state s — sq of the same equation. 

(c) Each state s is characterised by an its own probability density function (PDF) 
ips : — > M.'^ of "transition events": 

(1.5) ijjsit), with I ilJs{t)dt = 1. 

Jo 

(d) When a transition event occurs, the dynamics switches instantaneously from 
the state j (Aj) to a new state i (Ai), given randomly by the transition prob- 
ability matrix (or transition measure): 

(1.6) 

The position X{t) of the process is not affected when the state switches. 
Assumptions (|1.4p . ()1.5|) and ()1.6|) define the three local characteristics of the 
PDP. We see that Eq. ()1.4|) can be integrated as an ordinary differential equation, 
provided that not any switching event happen inside the integration interval. There- 
fore, with the exceptions of switching times, the process is deterministic, continuous 
and composed of pieces of solutions of Eq. (|1.4p . Anyway, the whole resulting process 
X{t) is not deterministic, it represents a random sample path on a probability space. 
1 

The statistical description of the process is given by the unknown functions 
Fs{x, y, t): each represents the probability to find the process X{t), in the state s, at 
time < in a position less than x, being past the time Y since the last switching event. 
Formally we write: 

Fs{x,y,t) := P{X{t) <x,y <Y <y + dy, state = s) 

where P is a probability measure of a probability space for the process. If we are 
interested only in the position x of the process, we can integrate over all values of y, 
and the distribution function for the process, regardless the time y, reads as: 

rt-to 

(1.7) J^s{x,t) := / dyFs{x,y,t). 

Jo 

With the further hypothesis X{t) £ ft, we have: 

(1.8) d^Fsix,y,t%^gn=0 

since there is a null probability for the process to be outside the interval f2. Besides, 
the probability measure have to be conserved during the evolution, so that: 



s 




For our pourposes, we do not need to specify what the abstract probability space is. 
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and 



(1.10) lim y^J^six,t) = 0, Wt 

s=l 

have to be satisfied. This three last equations are boundary conditions for p.ip . that 
complete the definition of the problem we approach to treat here. 

The function X{y), named hazard function (or hazard rate), is related to the 
statistics of the PDF switching times (|1.5p by: 



(1.11) A.(,) - 



It represents the probability per unit of time that a transition event will occur, i.e. a 
transition rate, having past the time y since the last event. The explicit dependence 
of A on ?/ makes both the statistics of the switching events and the process X{t) be 
non-Markovian, so that y plays the role of memory. 

The main aim of this article is to solve Eq. (|l.ip . jointly to all the above men- 
tioned boundary conditions, by a finite difference scheme of the first order and prove 
convergence of the numerical solution. We note that numerical methods for solving 
linear hyperbolic PDEs with non-constant coefficient, such as (jl.l[) . are well known 
in literature [H H] , but what makes this problem a little special is the boundary con- 
dition (jl.3p : the value of the unknowns Fg, on the boundary y = 0, depend on an 
integration of the Fs over the interior of the domain. This means that the numerical 
scheme for Eq. (|l.ip have to be supported by one for (|1.3p . so that conditions for 
convergence of the numerical solution have to be investigated again. As a result we 
found a Courant-Friedrichs-Lewy (CFL) condition for ensuring linear convergence. 

The secondary, but not of minor importance, aim of this article is to provide 
a connection bridge between PDPs as known by experts of the field and, as above 
mentioned, the same processes as known by others, who apply them to modelling in 
several areas of science and engineering. Here we give a sample of quotas, for which 
PDPs can be concerned by others, grouped in two categories: diffusive processes 
and systems having an equilibrium. We mention: anomalous diffusion |llj . reaction- 
diffusion [H], scattering of radiation [19], biological dispersal [27], for the former 
category, and non-Maxwellian equilibriums [HI HI [HI [IZl [HJ [JT] , diagnostic techniques 
for semiconductor lasers [2^ , filtered telegraph signals [28l [20] , harmonic oscillators 
[23] . ecological systems [22], for the latter. Many of the models involved in such 
references, concern the application of a two-state noise to a dynamical equation. 

The common end of all these researches, consists in extracting statistical prop- 
erties from process governed by that equation. Generally an approximation method 
can be applied to the original model: such as by the projector technique [51 [5]. 
by a "coarse-grain" technique, (see, e.g., [55]), an asymptotic analysis (see, e.g., [TT]). 
However not always these techniques provide a satisfactory description. In some cases 
an exact analytical result can obtained as in Refs. [28l|20l[25] and more recently by the 
characteristic functional method 12J. Obviously, computations can also be performed 
by Monte Carlo's simulations, but, at the best of our knowledge, few or nothing it 
has been devoted to a direct calculation of the time-dependent distribution function 
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comprensive of an explicit memory variable. Concerning this, we remark that the 
main alternative is based on the inclusion of supplementary variables [51 \TE\ , that 
turns PDP into Markovian, i.e. a memoryless process. 

In the next section we provide an example that emphatizes the connection between 
PDPs and models with dichotomic noise, and a conjecture that ensures the existence 
of a stationary solution of Eq. p.ip . In Sect. [3] we establish the numerical scheme. 
We introduce definitions in Sect. [4] and in Sect. [5] some theoretical results about 
the related convergence. In Sect. Hlwe present numerical results to two problems for 
which an analytical stationary solution is known in closed form, and verify the stated 
convergence properties. 

2. Explanatory example. Let us consider a dissipative process X{t) subject 
to a noised input ^{t), described by the equation: 



If ^(t) is taken as the random telegraph signal, Eq. (|2.ip act as filter, and X{t) is 
referred as filtered random telegraph process [2H1 [Hj ■ The same equation is elsewhere 
referred as Langevin equation [HI US] subject to a dichotomous noise, ^(i) alternately 
takes on values ±1, with an exponential (or Poisson) statistics for the transition events 
(jl.Sp : iP{t) = /ie"''"^, where is the average time (r) between transitions. The 
process X{t) results composed of pieces of increasing and decreasing exponentials. The 
statistical properties of the process X{t) can be found by the associated probability 
density distributions p^{x,t), governed by a Liouvillc-Mastcr Equation ^32. ,20. ,S\: 



Now let us to see the same process from the point of view of PDPs. The exponen- 
tial statistics for ip{t) makes the process of transitons be Markovian and the hazard 
function constant: X{t) = jjL. Eq. turns into: 

dtFs{x,y,t) + As{x)dxFs{x,y,t) + dyFs{x,y,t) ^ ~^Fs{x,y,t). 

By integrating this equation over all the values of y, we get: 

dtTs{x,t) + As{x)dxJ^s{x,t) - Fs{x,0,t) = -fiTs{x,t), 

having used the property Fs(x,y,t) = if y > 0, since the process is memoryless. 
From Ea. p.3p we have: 



(2.1) 



dX 

dt 



(2.2) 



{ 



dtP+ - {X- l)(9a;P+ = (1 - i-i)p'^ + i-ip 
dtp^ ~ (x + l)dxP^ = + (1 - 




and inserting it into the previous: 



s 



dtTs{x,t) + As{x)dxf^s{x,t) =^^(q^j- - Ssj)Tj{x,t) 



If S* = 2, with the known functions: 



(2.3) Ai{x)^{l-x), A2{x) ^ -{1 + x), 
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and transition measure: 

(2.4) qii = q22 = 0, qi2 = q2i = 1, 

provided that p^{x,t) — dxJ-i,2{x,t), we obtain just the equation (|2.2p . This show 
the connection between PDPs and processes driven by dichotomous noise. 

2.1. Remarks on equilibrium solutions. In what follows we focus our atten- 
tion on solutions F(x, y, t) having an equilibrium, but we presume that the numerical 
scheme can be extended to diffusion processes too. Conditions for the existence of 
equilibrium solution can be conjectured by using simple dynamical arguments |25(ll0j. 
If all dynamical equations (|1.4p own only attraction points and all these are contained 
into the intersection of the basin of attraction of each, then a process starting from this 
region will never escape. Whence, there should exists a region SI where the process is 
confined and a stationary distribution J-c,q{x) = limj^^oo -^(2;, t) exists. 

3. The finite-difference scheme. In this section we show the numerical scheme 
to solve Eqs. and (|1.3p based on a finite difference method of first order. 

For the shake of simplicity we take to = and the domain of Fs becomes T) := 
{fl X [0,T] X [0, r]}. It is convenient to perform the numerical integration along 
the characteristic lines ^ = t ~ y. With this new variable, we define the unknowns 
4>s{x,y,S,) = 4>s{x,y,t ~ y) :— Fs{x,y,i), so that Eq. p.f p transforms as: 

(3.1) Ai,{x) dx(j)s{x,y,£,) + dy(j)s{x,y.X) = -^s{y)(l}s{x,y, £,)■ 

This equation is valid for < ^ < t and < y < t. The initial condition is given on 



(3.2) M^,y,Ok=-y 

With the new variable we get (j)s{x, 0, £,)\i=t 
Eq. (|1.3p becomes: 



= Fo.s{x)Siy) 

= Fs{x, 0, t), and the boundary condition 



(3.3) (j)s{x,Q,t) ^^qsi I dy (j){x,y,C)\e,=t~yMy) 

1=1 -^0 

We will assume that similar conditions of Eqs. (|I.7p . (|I.9p and (|I.IOp are satisfied for 
(/>s(x,y,^), and a stationary solution exists. 

On the domain I?, we introduce a uniform mesh: 

^^■4) I ,-„^o,...,7V, N = T/At 

with step size Ax and Ay — At, so that we define the discrete known functions as 
''Ak '■— Ai{xk) and 'Aj :~ Xi{yj), and the discrete solution: 

'Ffe"^-, n = 0,...,iV, j<n 

as an approximation of Fi{xk, yj, tn) at the mesh points. 

The change of variable i, ~ t — y corresponds to the following discrete mapping 
on the mesh: 



(3.5) 



(fc, j,n) = [k,j,i)\i=n~j, 
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therefore we get the following relation: 

(3.6) Fi{xk,yj,tn) = (t)i{xk,y„t,,-yj) = (t>i{xk,y,,£.i) ^ ^F^^ = ^^7^ = % 

between the discrete solutions. Here the index i identifies the set of mesh points lying 
on the characteristics lines. 

The numerical scheme is obtained by discretizating both equations (|3.ip and (|3.3p . 
We apply upwind to the first, and get: 

(3.7) = % ~ '^.^ ( - ''/'U.-i,,) - % 

i^O,...,N 

where v — 1 ii < 0, and i' ~ Q ii ^A^ > 0. The boundary condition (jl.Sp is 
included by requiring that '(/ig^- = and '^tv^-i j — ''^a'^ j- 

For the second, we substitute integral with a quadrature scheme: 

(3.8) %o = Ay 5] qsi E % ^ > 

1=1 j=0 

where Wj > is a sequence of weights. 

The integration proceeds as follows. Given the initial condition '0° o ~ i^f^ i 0^ 0) 
= i^/(a;fe, 0, 0), Eq. (13. 7p allow us to calculate all the '0^^- starting from j = 1 up to 
N , for the fixed characteristic line i — 0. In general given the values on the boundary 
j = 0: '0^Q, we can find all ''(f>l.j starting from j = 1 up to — i, for a fixed 
characteristic line i (see curved arrows of Fig. dSj])). But starting values '(^^^^ for 
upwind are unknown and have to be estimated by using the boundary integration 
p.Sp (see vertical dot-dashed arrows of Fig. (13. ip ). This is a system of equations for 
the unknowns ''ipl, q . When all '0^^ are known, the discrete distribution function can 
be retrieved by 'i^^^ = '■6lj\i-„-j. and also the discrete countrepart of (|1.7p : 

n 

(3.9) '.f," :=^«f^'i^,",. Ay, 

can be estimated by a quadrature formula of weights wj"^ . 
4. Preliminary definitions. 

4.1. Global errors and convergence. We are interested in how well '0^^ , 
^F'^j and ''J-J} approximate the corresponding analytical solutions. We consider the 
global pointwise error: 

(4.1) := Mxk,yJ,^^) - % 

for the transformed solution. The same value defines the global error for ''F^j under 
the discrete mapping (|3.5p . 

In order to prove convergence, we introduce norms for measuring errors. For 
spatial Xk and memory variable yj the cx)— norm is used. The discrete 1-norm for 
the states s of the system is the natural choice, because of the conservation of the 
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Fig. 3.1. Representation of the integration scheme on the regular mesh of the {t,y) plane. Full 
curved arrows: upwind step of Eq. {S.Tjj . Dot-dashed vertical arrows: quadrature of Eq. (g.gj) 



probability of Eq. (|1.9p . For convenience of notation we define the global error for 
the state / at time ti and time memory j/j as: 

(4.2) ^ :=max| ^,,1 
and the global error regardless states as: 

s 

(4.3) ll^^lli:=E 

1=1 

We say that ^FJ^j converges to Fi{xk,yj,tn) in the norm || • || if: 

(4.4) Pi = max UK 111 ^ 0, as Ax, Ay -^0. 

j 

The global error for the distribution function (|3.9p is defined as: 

(4.5) 4" := H^k,U,) -TJ:= dyY,Fi{xk,y,U,) - Ay^z;f) ^ 'F^,^ 

•J" I j=0 I 

and the associated convergence is stated by: 

(4.6) ||^"Hoo = max|i?^'| 0, as Ax, Ay ^0. 

k 

4.2. Local truncation error and quadrature error. As usual [HIS], the local 
truncation error is defined by inserting the true solution (j)i(x,y,^) into the discrete 
scheme of Eq. (13. 7p . i.e.: 

I t ^ (t)i{xk,yt + Ay,^0 - (l)i{xk,yi,^i) 
Ay 

(j)i{xk + i^Ax,yi,£,,)-(j>i{xk + {i^-l)Ax,yi,£,,) 
(4-7) -Mxk) 

-'k{yj)(t)i{xk,yi,Ci)- 



By evaluating the rest with respect to '0^^-, we get: 

(4.8) = ^Ax {\AiixkMM7jk,yj,Cd - adld^i{xk.rij,S,,)) 

where a := Ay/Ax^ and rj^, rjj are unknown points. 

The quadrature error committed from (|3.8p for the evaluation of the boundary 
integral ()3.3p is defined as: 

(4.9) % := M^UAy^^^{M^k,f|J,^^) MVi)) 

where Mi are some constants that can depend on i, and fjj are unknown points of the 
local integration interval. /3 defines the order of repeated quadrature formulas (|3.8p . 
The quadrature error committed from (13. 9p for the evaluation of p.7|) is defined 

as: 

(4.10) := MntnAthlFi{xk,Vj,t„), 

where, as for the previous error, Mn are some values that can depend on n, and fij 
are unknown points. f3 defines the order of (|3.9p . 

5. Analysis of convergence. In this section we show first a lemma for conver- 
gence of the numerical solution for the transformed equation (|3.7p . then proof 
a theorem for the convergence order of the numerical solution J-'J} to the distribution 
function !F{xk,tn)- Proofs are based on classical arguments by finding bounds for 
global errors. 

Lemma 5.1. Let 4'i{x,y,£,) e C^''''^(I?) be a solution of Eq. i3.1\) under the 
boundary conditions !i3.2\} . with (3 = max{/3,2}, max; ||Ai(a;)||oo ^ M, and 
max; ||Ai(j/)||oo < Lu, for x € n and I € {1, . . . ,S}. Let {xkTyj,tn) be an uniform 
mesh on "D defined in {3.4^ , of step sizes Ax and Ay = At. Let ^(f)\j be the numerical 
solution resulting from the scheme as defined in Eqs. {3.1^ and fi3.8\) . under the 
transformed mesh of Eq. 

If the Courant-Friedrichs-Lewy (CFL) condition 

(5.1) ^^<(£+^")" 
is satisfied then: 

1. Given the error ||i?olli ^.t boundary y — 0, the error computed at time 
.step ti+rn along the characteristic is bounded by: 

(5.2) < ||i?^||i(l - LuAyy- + mAy\\£% 

where \\£'\\i =0{Ax), Vi. 

2. If Ay^^ > max,(w^'') max; Ai(0), there exist constants Ld, K , such that the 
error computed at time ti along the boundary y — is bounded by: 

(5.3) < p||iexp(i^t,) + 11-^11 i-^^exp((X-Lrf)t,) 

where \\R\\i = 0{AxP) and \\8\\i := max, ||£^||i = 0{Ax). 
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This lemma states that the numerical solution '0^^ converges to analytical with 
a linear order. Errors calculated for (j) are same for F . Being interested in finding the 
probability regardless of memory and state at fixed time t„ , we search an estimate for 
the error H-E'^ljoo as defined in Eq. (|4.6p . 

Theorem 5.2. Let the hypothesis of the Lemma \ 5.1\ he satisfied, then for the 
problem of Eq. hLl\) the numerical solution k3.9\) . obtained as discussed in Sec. 
converges to analytical Jj. ?| ) in the sense with order ©(Ax). 

Proof. We substitute the integration of Eq. (|1.7p over the memory state with a 
sequence of quadrature formulas of weights wj""* : 



(5.4) / dy^i^K2;fc,y,t„) = Ay^«f '^FKxfc,y„i„)+^ 'i?^' 

•'^ I j=0 I I 

so that: 

j I I j I 

Let w*-"-* := maxj and R" := J2i max^ | from the triangle inequality 

II ^" lU < Ay E E 1^' (^'^ ,yj,tn)- 'Fl^ I + i?" 
by using the discrete mapping (|3.5p and the definition of the error (|4.3p . we get: 



|^"||oo <«("^AyEll^^'lll + ^"■ 
i=o 



By inserting Eq. (j5.2p we find: 



|^"||oo < «(")AyE(ll^o"'lli(l - L«Ay)^' + jAy|l£ii) 

3=0 



and 



n _ f2 

P"||oo <^^^"^AyE ll^o"lli(l-inAy)"-"+z;(")||£||i^ + 7?" 

m=0 

where ||f ||i := maxj ||£'"~^||i. This inequality relates the searched probability error 
to errors along the boundary j — Q. Now we insert the second result of the previous 
lemma stated by Eq. (|5.3p and find the order of the error for vanishing Ay. From the 
summation we get: 

" / f2 

E(l - Lu/^VT-^ Pllie^*- + m,K-^e^^-^-^'^ 

that is of order: 

Pill , ' , , +l|g||li^ 



{K + Lu)Ay " 2iK + Lu~Ld)Ay 
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Finally, we get an order of convergence for the error: 



(5.5) p"|U < ^(")p|U "^P^,^^"^ ^ + i?" 



□ 



6. Computational results. In what foUows we present results of the numerical 
scheme, apphed to the examples considered in Ref. [53]. For quadrature of Eqs. 
and (|5.4p . we adopt the rectangle scheme: 

(6 1) w^"^ = v^"^ ^ r for J = 

> ^1 1 1 1 for j > ^ 



whose quadrature error (|4.9p is: 

'^i = ^i»Ay9j,(0,(xfc,%,,ei)A;(?7j)). 

Note that despite the low order {(3 = 1) of approximation for quadrature, the global 
error of the numerical solution is not degradated, because the same order apply to 
upwind. The choice = makes the equation (|3.8p be explicit. This can be 
consistent, because for vanishing Ay the contribution to integral is also vanishing for 
a limited ^4>\ q . For the explicit scheme the computational cost can be evaluated as 
follows: at time step i all upwinds take 2NkSi operations, the boundary quadrature 
takes NkS'^i. By summing over i we get: 

Computational Cost « NkS'^N'^. 



The Cauchy problem for starting the numerical integration is set, according to 
Eq. (fr2|) . as follows: 

r 0, fc<0 
'F^o^i 0.5/(5Ay), k = 
(6.2) [ l/{SAy), k>0 

^FO. =0, j>0 

for all s = 1, . . . , S". This choice is an approximation of (|1.2p with Fo,s{x) = H{x), 
where H{x) is the Heaviside function. Such Cauchy conditions for the Liouville- 
Master Equation correspond to having placed the process X{t) at the initial position 
Xq = 0, to an equiprobable random initial state , having spent the time zero (i.e. S{y) 
in (|1.2p ). This is a common choice when studying the motion of a particle subject to 
a random fluctuating force, but is not a good mathematical hypothesis for applying 
Lemma [5.11 However, it is well known the upwind methods tends to regularise the 
solution (numerical viscosity) around discontinuities [SllTj, and, for the problems we 
are approaching to solve, a unique stationary solution exists regardless the initial state 
of the process. In this way we are enabled to use such "non-regular" Cauchy condition 
for our numerical convergence tests. 'F^^ are then integrated by (|5.4p . 
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Fig. 6.1. Six snapshots of the temporal evolution of p^^ for the Langevin equation driven by 
Poisson distribution time intervals of Sec. \6.1\ (horizontal = 0,4,8, 12, 16,20;. 

We are interested in plotting the density probability function of the process, 
regardless memory and states, defined as: 



(6.3) p{x,t) ■.^d^T{x,t) := d^^Ti{x,t). 

1=1 

The discrete version of this operation is: 

(6.4) p: := - ^fe ) := 7^ E( '^^+1 ^ 



that is the first order right-derivative of the numerical distribution function. 

6.1. RC-filter subject to Markovian process (Poisson PDF). In Fig. 
are plotted six snapshots of the temporal evolution related to the four states process: 
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Fig. 6.2. Six snapshots of the temporal evolution of p^ for filtering of dichotomous noise driven 
by the McFadden distribution time intervals (horizontal axis: x^). 



As{x) = -^sX + Ws, studied in [71. Parameters are: Wq = 1,Wi^ -1, W2 = 2, W3 = 
—2, Xg — 0.2 and 7^ = 0.1 Vs. Results are comparable with that of the cited 
reference. Note the broadening of the four peaks due to the numerical viscosity of the 
upwind method. 



6.2. Filtering of non-Markov dichotomous noise with McFadden inter- 
val PDF. For this example the process is described by the Langevin equation (|2.ip . 
with functions of Eq. (|2.3p and transition measure of Eq. (I2.4p . Intervals between 
switching time have the McFadden PDF: ^p{t) — 3e^*(l — e^*)^. The equilibrium 
density distribution has the form [5S]: 

3 

(6.5) peq{x)^ \imp{x,t)^—{7 + x^), \x\ < 1 

t^OQ 44 

12 



and its integral: 



(6.6) Ten{x) = \im T{x,t) ^ ^{x^ + 21x + 22), \x\ < 1 

The hazard function related to the density for switching intervals is: 

3(l-e-*)2 



m = 



(1 



We note that is A(0) — 0, so that the error committed from the choice Wq'^ — (see 



(16. is further improved. Beside it is X{t) < 4/9 and the convergence lemma give 
us more guarantee that the errors does not grows fastly. 

Here the grid step sizes are Ax = 0.002 and Ay = 8 • 10^^. Integration starts 
with a concentrated initial density (|6.2p and stops at time T = 7.37. 

In Fig. 16.21 we see six snapshots of the numerical solution of the PDF (|6.4p . 
At time t — the density of the process is concentrated to a; = 0, then two peaks, 
corresponding to the two dynamical states, moves towards the attraction points x±l, 
and at last the stationary distribution appears. 

6.3. Filtering of non-Markov dichotomous noise with "gamma" interval 

PDF. For this example the process is described by the same Langevin equation as 
that the previous one, but the intervals between switching times of (,{t) of Eq. (|2.ip 
are taken as the gamma density ^p{t) — /i^te"'^* for both states [55]. Provided that 
/i = 1/2, the equilibrium solution for the total density distribution function is: 

(6.7) pe,{x) = lunp{x,t) = Ml_^(i_^2)-i/2 1.2;2)^ 

= {l + ^/^)/A, \x\<l 



r 



in which 2F1 is a hypergeometric function. 

The hazard function (jl.lip related to iplt) is: 

Mt)- '''' 



l + ^J,t 



We see from the convergence theorem that the error does not grow so fastly, because 
the maximum value of A(i) is Xmax = m/c We perform the numerical integration on 
a mesh with spatial discretization step Aa; = 0.004 and temporal step At = 0.0015. 
Integration starts with a ()6.2p and stops time T = 12, where the equilibrium is 
supposed to be reached in good approximation. 

In Fig. 16.31 are plotted six snapshots of the total density distribution (|6.4p . The 
evolution behaves as the previous example, with the exception of singularities at a;± 1, 
for the equilibrium. 

6.4. Convergence tests. Since for the above mentioned problems we know 
two analytical results, we can calculate the global error ||-E||oo for the stationary 
distribution functions of Eqs. (|6.6P and (|6.7p . The analytical solution of Eq. (|6.7p 
is evaluated by using MATLAB(gO with libraries for calculating the Hypergeometric 
function [30]. We calculate the solution with the numerical scheme until the time 
T = 20. At this time we experienced that the stationary solution is reached. This 



^MATLAB® is a trade mark of "The Matworks, Inc." 
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Fig. 6.3. Six snapshots of the temporal evolution of p^ for filtering of dichotomous noise driven 
by the gamma distribution time intervals (horizontal axis: Xk). 



integration is repeated for some spatial step size Ax, with the temporal step size 
constraint a = Ay/ Ax — 0.9 (M + L„Aa;)~^, satisfying the CFL condition (|5.1|) . 
In Fig. ((Ol) we show the global error Ek plotted for Ax = {0.1,0.04,0.008} for 
the McFadden intervals. We see clearly that the maximum error decreases as Ax 
decreases. In Fig. (|6.5p we show the same test for gamma distributed intervals. Also 
here the error decreases, but it shows a sort of divergence near x — ±1. 

In order to stress convergence, in Fig. (|6.6p . we plot the error ||i?||oo versus the 
step size Ax. We see that the McFadden's data's have unitary slope, i.e. linear con- 
vergence, as we expected from Theorem (|5.2I) . Instead gamma's data's are arranged 
with approximately 1/2 slope. This can be explained as follows. We know [29] that 
the PDF of Eq. ((O)) is U-shaped having (1 - x^)~^/^ ln(l - x^) infinities near x ± 1. 
As (5x := 1 — X approaches to zero, the second spatial derivative of J^e„ behaves as: 



d^Peq = dlTe, oc + ln(2<5x)). 
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Fig. 6.4. Global error Ei^ between the exact solution of Eq. 16. 61) (McFadden PDF) and 
numerical calculated for Ax = 0.1, Ax = 0.04, Ax = 0.008. 
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Fig. 6.5. Global error Ef^ between the integral of the exact solution of Eq. \6. 7] ) (gamma PDF) 
and numerical calculated for Ax = 0.1, Ax = 0.04, Ax = 0.008. 



Being Ai{5x) ~ Sx and by considering Ax — 5x, we get for the error: 



\E\\ 



ln(2 Ax)), 



that explain what we see from the results. If we remove the endpoints from the 
measure of the error, we see from "gamma 90%" of Fig. 16.61 that, e.g., for the 
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Fig. 6.6. Test for global GTTOT ||£^||cxD ^S. TTlGsh Step siZG /\x , 



interval x £ [—0.9, 0.9] the linear convergence order is recovered. 

Appendix A. Proof of Lemma 15.11 

Proof. 

We study convergence in two parts: the first concerns the accumulation of the 
error along the characteristic lines of Eq. (j3.7p : the second concern the error cumulated 
on the boundary integration of Eq. (|3.8p . 



Let 'e^ be the global error (|4.ip between exact and discrete solutions. By 



considering the standard procedure O [4] that puts in relation local and global errors, 
we find from Eqs. ([XT]) and 

(A.l) = '4,, - 'A,a{ - H+.-i.,) - % H,,A2/+ H,Ay 



where 'e^^- is defined in ()4.7p . Wc start our analysis from Eq. (jA.ip . fixed i and given 
j we find 'e^. j^^^. The error for 'e^. g is found from the boundary integral condition. 

The first step is to limit the error along the characteristic line i. From Eq. (jA.ip 
we study the case > {i^ — 0): 



Li 



< |1 - 'Aka - 'A, Ay I | + H-iJ + I '4, 1 Ay 



Let ^Ej = maxfe | 'e^ j \ and '£j = max^ | 'e^^ | then: 



If we set 1~ 



< (|1 - 'Aka- %Ay\ + 'Aka) 'E] + fe]Ay 
^\j/S.y > VZ,j, fc, i.e. the CFL condition, we get: 



'4,,+i I < (1 - %^y) 'E} + '^j Ay Vfc, J, I 
16 



because of it is valid for all k, we can write: 

(A.2) 'E]^, < (1 - 'A, Ay) 'E} + fejAy 

that is verified for j = 0, . . . , n and i = 0, . . . , n. Now we have to find ^Eq^^. 
The case v = 1 gives the same result: 

I 'eL,+i|<|l+ 'Aka- %Ay\\ 'dj - H+iJ + I '-^jlAy 

and for the maximum error on a;: 

I <{\l-\'A,\a- %Ay\ + \ 'A,\a) 'E] + fejAy 

Let 1 — I ^Ak\a — ''XjAy > 0, we get Eq. (jA.2[) . so that it is valid independently by 
the sign of Ai{x). 

By iterating this expression we find: 

(A.3) < n™ o(l - 'a, Ay) 'Ei, + {rn + l)Ay fe' 

where ' ~ max., 'fj. When the norm over all states is considered: 
(A.4) < ||i?^||i(f-i„Ay)™ + mAy||r||i 

Being 'Aj > we have good chances to get an upper bound to the error for 
increasing time. 

Now we study the second step, i.e. the error 'e^ q (or ''Eq) along the boundary 
condition. Here j = then n — i. 

From Eq. (|3.3p as calculated with quadrature on the exact solution y, ^), we 
can write: 

S I i 

1=1 \3=0 

Subtracting side by side this equation and (|3.8|) : 

s 



1=1 \]=n 



and by applying the triangle inequality: 

I V,,ol < Eg.; I E Ay«;« 'A,| H"'! + | % 
1=1 \j=i 

Let = maxfc | 'i?^ and being qsi > 0, wj*"* > 0, we get: 

(A.5) ^E^o < E E ^y^f + 

1 = 1 \3=0 

Let 'A = miuj 'Aj, from Eq. (|A.3p is: 

'E]^, < (1 - 'AAy)^+i % + {] + I) Ay 
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and inserting it in (jA.Sp . we find: 



(A.6) ^E^^ <^g,J Ay^zwf 'A, [(1- 'XAyY + jAy + 'i?' 

1=1 \ 3=0 

This inequality puts in relation the error on the boundary condition {j = 0) with 
those at early times. 

Let w)*-*^ = inaxj>i and 'A — inaxj>i 'Aj, then: 



^i?o <J2l^^ {^y^^^^ '^E[(l- 'AAyy'Eir'+jAy'r-^] 

1=1 \ 3=1 

+Ayw^;^ X % + 'R' 

valid for i > 1. We have introduced the sum starting from j = 0, so that, in order to 
find ^i?o , we have to solve a system of equation. However for convergence we do not 
need to solve it at all. By moving to the left hand side and summing over the 
states, we find: 



^ ^El^Y.'i^iAyw^^^X'ElK 



s s,l 

S / i 

^^q,, Ayu;« 'A^[(1 - 'XAy)' + jAy 'f'"^"] + 'ff 

s 1 = 1 \ j=l 

By using the fundamental property of stochastic matrix: 
^(1 - woLoAy) % < 

i 

3 = 1 I I 

where wq — max^ Wq \ Lq = max; Aq. Let a := 1 — woLoAy > 0, we get: 

I 3=1 I I 

where < a < 1. Let i„ = max; 'A, Ld = min; 'A, = (1 — LdAy) provided that 
from CFL it is surely LdAy < 1, we have: 

i 

(A.7) \\El\\, < a-^Ayw^^ uY.^V\\Er\\, + jAy\\£'-^\,] + \\R%, 

and with the summation index change m = i — j , it becomes: 

6-'||£;^||i < KAy ( Ay!(l±ll||£-||i + ^ ) + 

\ m=0 / 
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where ||5||i = max^ ||£'||i, and K = maxi(a -'^w^*' L„). From the discrete Gronwall 
lemma (cf. in Th. 1.5.4 e Corol. 1.5.1 of [1]), we obtain: 

b''\\E',\\i < (^KAy^'-^^m, + (1 + KAyy. 

Let ||-R||i = maxj ||i?''||i, we find: 

< \\RUl + KAyy + \\EhKAy'^-^^[b{l + KAy)Y 

and finally: 

(A.8) WE'oh < Pill exp {KU) + miiK*^ exp {{K - L^) U) 

□ 
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